#include <iostream>
#include <ctime>
#include <x86intrin.h>
//#include
using namespace std; 

double compute_pi_naive(size_t dt){
	double pi = 0.0;
	double delta = 1.0/dt;
	for (size_t i =0;i<dt;i++){
		double x = (double)i/dt;
		pi += delta /(1+x*x);
	}
	return pi*4.0;
}

double compute_pi_omp_avx(size_t dt){
	double pi = 0.0;
	double delta = 1.0/dt;
	__m128d ymm0,ymm1,ymm2,ymm3,ymm4;
	ymm0 = _mm_set1_pd(1.0);
	ymm1 = _mm_set1_pd(delta);
	ymm2 = _mm_set_pd(delta,0.0);
	ymm4 = _mm_setzero_pd();
	for (int i = 0;i<dt-2;i+=2){
		ymm3 = _mm_set1_pd(i*delta);
		ymm3 = _mm_add_pd(ymm3,ymm2);
		ymm3 = _mm_mul_pd(ymm3,ymm3);
		ymm3 = _mm_add_pd(ymm0,ymm3);
		ymm3 = _mm_div_pd(ymm1,ymm3);
		ymm4 = _mm_add_pd(ymm4,ymm3);
		
	}
/*
	double tmp[2] __attribute__((aligned(16)));
	_mm_store_pd(tmp,ymm4);
	pi += tmp[0]+tmp[1];
*/
	ymm4 = _mm_hadd_pd(ymm4,ymm4);
	_mm_store_sd(&pi,ymm4);
	return pi*4.0;
}
/* 

double compute_pi_omp_avx(size_t dt){
	double pi = 0.0;
	double delta = 1.0/dt;
	__m256d ymm0,ymm1,ymm2,ymm3,ymm4;
	ymm0 = _mm256_set1_pd(1.0);
	ymm1 = _mm256_set1_pd(delta);
	ymm2 = _mm256_set_pd(delta*3,delta*2,delta,0.0);
	ymm4 = _mm256_setzero_pd();
	for (int i = 0;i<dt-4;i+=4){
		ymm3 = _mm256_set1_pd(i*delta);
		ymm3 = _mm256_add_pd(ymm3,ymm2);
		ymm3 = _mm256_mul_pd(ymm3,ymm3);
		ymm3 = _mm256_add_pd(ymm0,ymm3);
		ymm3 = _mm256_div_pd(ymm1,ymm3);
		ymm4 = _mm256_add_pd(ymm4,ymm3);
		
	}
	double tmp[4] __attribute__((aligned(32)));
	_mm256_store_pd(tmp,ymm4);
	pi += tmp[0]+tmp[1]+tmp[2]+tmp[3];
	return pi*4.0;
}
*/ 
int main(){
	clock_t start,end;
	size_t dt = 128000000;
	double result1,result2;
	
	start = clock();
	result1 = compute_pi_naive(dt);
	end  = clock();
	cout<< result1 <<endl<<end- start <<endl;

	start = clock();
	result2 = compute_pi_omp_avx(dt);
	end  = clock();
	cout<< result2 <<endl<<end- start <<endl;

}
